Geometrical defects in two-dimensional melting of many-particle Yukawa systems 
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We perform Langevin dynamics simulations and use polygon construction method to investigate 
two-dimensional (2D) melting and freezing transitions in many-particle Yukawa systems. 2D melting 
transitions can be characterized as proliferation of geometrical defects — non-triangular polygons, 
obtained by removing unusually long bonds in the triangulation of particle positions. A 2D liquid 
is characterized by the temperature-independent number of quadrilaterals and linearly increasing 
number of pentagons. We analyse specific types of vertices, classified by the type and distribution 
of polygons surrounding them, and determine temperature dependencies of their concentrations. 
Critical points in a solid-liquid transition are followed by the peaks in the abundances of certain 
types of vertices. 
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I. INTRODUCTION 

For a few decades melting and freezing transitions 
in two-dimensional many-particle systems have been in- 
vestigated in a variety of experimental and compu- 
tational studies, without reaching a definite conclu- 
sion regarding its nature. According to the most 
widely accepted Kostcrlitz-Thouless-Halperin-Nelson- 
Young (KTHNY) theory, melting of a single two- 
dimensional (2D) crystal occurs via two continuous phase 
transitions, first from a solid to hexatic phase, and then 
from the hexatic to an isotropic fluid [TJ [5] . The theory 
also predicts, what it is proliferation of unbound topo- 
logical defects — dislocations and then free disclinations 
— what plays a crucial role in 2D phase transitions and 
breaks positional and orientational order. Some exper- 
iments and theoretical studies, however, suggest a first- 
order grain-boundary-induced melting scenario in poly- 
crystalline systems [3l 0] . 

Although in 2D a true crystalline order can not sur- 
vive at finite temperatures kT > [5J |B], a quasi-long 
range translational order is observed at the conditions of 
strong coupling, for example, in complex plasma layers 
[7] or charged colloidal suspensions [8\. Over the years 
a broad range of empirical criteria was developed to ac- 
curately determine melting and crystallization points [9] . 
Perhaps one of the most famous examples is the Lindc- 
mann criterion, which has been applied extensively in 
three-dimensional melting and freezing, while its gen- 
eralized version has been used in some studies of two- 
dimensional transitions [10 . Other methods frequently 
make use of topological defect fractions, bond orienta- 
tional order parameters, orientational and positional cor- 
relation functions as well as Einstein frequencies [11) . In 
a recent work, a polygon construction method by Glaser 
and Clark [T2TH^] was employed to characterize tran- 
sitions in a rapidly heated and cooled two-dimensional 
complex plasma experiment [15] . 

Systems of strongly correlated particles in complex 



plasmas are of particularly high importance in the ex- 
perimental studies of phase transitions. Complex plasma 
usually consists of polymer microparticles immersed in a 
weakly ionized gas, where distinct dust grains are known 
to interact through the Yukawa (Debye-Hiickel) interac- 
tion [TB] . Convenient time and length scales of these sys- 
tems allow for the direct optical observation of collective 
many-particle phenomena as well as accurate measure- 
ments of individual particle positions by the means of 
video microscopy [i~7Hl9] . 

In a recent experiment, the phenomenon of superheat- 
ing was observed in a solid state of two-dimensional com- 
plex plasma [20) . It was demonstrated, that during a 
rapid heating process, the concentration of defects can 
stay low and complex plasma can retain the properties 
of a solid even at the temperatures above a melting point. 
The same experimental results were later analysed by the 
method of geometrical defects [TS], originally developed 
by Glaser and Clark, which we use in the current work. 

Geometrical defects in a polygon construction method 
are identified by removing unusually long bonds in the 
triangulation map of particle positions, so that result- 
ing polygons have three or more sides. It was shown, 
that this method provides great sensitivity and unveils 
some interesting features, undetectable by the conven- 
tional analysis of topological defects [T5]. As another 
measure of disorder, the abundance of different kinds of 
vertices, grouped according to the type and order of the 
adjacent polygons, was suggested in the same work. 

In the present contribution we report numerical studies 
of two-dimensional melting and crystallization in strongly 
coupled Yukawa systems. Langevin dynamics simula- 
tions are performed to simulate gradual heating and cool- 
ing. Critical points are determined employing orienta- 
tional order parameters and topological defect fractions. 
However, the main motivation behind the current work is 
to present the method of geometrical defects and vertex 
fractions in the polygon construction as a sensitive tool 
to analyse the order-disorder transitions and characterize 
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the state of a 2D system. Furthermore, our findings re- 
solve the issue of the prominent peaks in the temperature- 
dependencies of certain types of vertex concentrations 
[T5] , by showing, that peaks correspond to the critical 
points in the initial and final stages of the order-disorder 
phase transition. 

In the following section we briefly describe the model 
system and simulation methods, as well as essential tools 
used in the analysis of phase transitions. Main results of 
the simulations and numerical analysis are presented in 
Section [TTT[ while Section |TV| summarizes the article. 



the friction coefficient /i is related to the Gaussian noise 
fj(t) by the fluctuation-dissipation theorem [22] 

(fi(t)W)) = ^kT Icl S ZJ S(t - 0, (4) 

where i, j € {1, . . . , N} and fcT re f is the desired target 
temperature in the units of e. In our simulations we use 
/i = 0.2. The Langevin equation is integrated numerically 
employing an impulse method of integration 23 . 

B. Analysis 



II. SIMULATION 

A. Model system 

A widely used approximation to describe interactions 
between particles in complex plasmas is the Yukawa 
inter-particle potential [16] 

Vij = Q 2 (47re r lJ )- 1 exp(-r u /A D ). (1) 

Here Q is the charge of a particle, is the distance 
between particles i and j, Ad stands for the Debye length, 
which accounts for the screening of the interaction by 
other plasma species. 

Strongly coupled many-particle systems with Yukawa 
interactions are fully characterized by two dimension- 
less parameters, namely the coupling strength T = 
Q 2 / {Airs^bk^T) and screening parameter k = &/Ad, 
where b is the two-dimensional Wigner-Seitz radius |21j . 
In the simulations presented here we set the screening 
strength to the constant value of k = 2. 

As a scale of length in our numerical simulations it is 
convenient to choose the Wigner-Seitz radius b, which is 
directly related to the areal number concentration of par- 
ticles, n = l/(7r£> 2 ). Therefore, the corresponding scale 
of energy is e — Q 2 / (Ane^b) and time is scaled according 
to the value of an inverse plasma frequency 



The model system consists of N — 2430 identical par- 
ticles in a 2D rectangular box of area 85.71 x 89.07, in- 
teracting via the Yukawa potential. Periodic boundary 
conditions are applied and, since the inter-particle poten- 
tial is short-ranged, the cut-off distance is set to r c = 8. 
Only particle pairs separated by less than r c are taken 
into account in the force calculation. 

We study order-disorder transitions in the model sys- 
tem by performing Langevin dynamics simulations with 
slow changes of temperature. Particle positions are up- 
dated according to the dimensionless Langevin equation 

¥i = -ViV(ri,...,T N )-(iri + fi, (3) 

where represents a randomly fluctuating Brownian 
force. In a thermodynamic equilibrium (fj(t)) = 0, while 



A common way of analysing the structure of 2D many- 
particle systems is calculation of a Delaunay triangula- 
tion, which yields a network of bonds connecting each 
particle with its nearest neighbours. A coordination 
number can be assigned to each particle, which is a num- 
ber of the triangulation bonds between the particle and 
its closest neighbours. The coordination number of a 
particle in a perfect hexagonal lattice is always equal to 
six. Topological defects are identified as particles with a 
different coordination number, usually five or seven. 

Two most common defect types are the disclination (a 
single particle with a non-sixfold coordination) and the 
dislocation (two connected particles with five and seven 
closest neighbours) [23]. Quite frequently defects orga- 
nize themselves in lengthy chains or grain boundaries, 
indicating a polycrystalline structure of the system. As 
an alternative, a Voronoi construction is sometimes used 
in the context of 2D dusty plasmas, where defects are 
identified as non-six-sided polygons [3J |2"5"H2"T] . 

To quantify the abundance of topological defects, we 
use the defect fraction DF, which is defined as a number 
of vertices with a coordination number other than six, 
normalized to the total number of particles N [2"7] . 

The polygon construction is a different way of charac- 
terizing defects in 2D systems O [13l [15] and helps to 
identify empty volumes in 2D liquids [25]. The authors 
of [T2] analysed bond-angle and bond-length probability 
distributions in dense 2D liquids. It was shown, that 
disordered regions of a liquid exhibit multiple peaks in 
a bond-angle probability distribution, with extra peaks 
corresponding to the square lattice. At the same time, 
triangularly-oriented clusters had a single peak and well 
expressed triangular order. The structure of a two- 
dimensional system therefore was described as a square- 
triangular tiling containing numerous tiling faults. 

A triangle, which is the only kind of polygon in the ini- 
tial triangulation of particle positions, is considered as a 
non-defective entity. To identify geometrical defects, cer- 
tain bonds are removed from the triangulation map, so 
that two polygons sharing a common bond are merged 
into one. Two possible approaches for the selection of 
bonds were suggested: either use a bond-length thresh- 
old, or identify a bond that is opposite to the unusually 
large angle between a pair of adjoining bonds. In our 
analysis we follow the authors of [HI [TS] and use the 
critical bond-angle of a — 75°. 
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The construction of a polygon map is illustrated in 
Figure [I] The first panel (a) shows the 2D triangula- 
tion map of particle positions, where most of the parti- 
cles are connected with six closest neighbours. Some ver- 
tices, however, have five or seven bonds and are marked 
by small triangles and squares. These particles are con- 
sidered as topological defects and are all the part of a 
lengthy grain boundary. Triangulation bonds marked by 
bold lines are facing bond-angles larger than the criti- 
cal value of a = 75° and therefore are selected for re- 
moval. The resultant polygon map is presented in the 
second panel (b). Geometrical defects are identified as 
non-triangular polygons, that is, quadrilaterals and pen- 
tagons. Although the defects appear in a close vicinity 
of the grain boundary, there is no one-to-one correspon- 
dence between the polygons and topological defects. 




and the abundance of geometrical defects we use four dis- 
tinct order parameters P n = N n /(2N) (n = 3, 4, 5, 6), 
defined as the number of polygons with n sides N n , nor- 
malized to the doubled number of particles 2N. 

Figure [2] provides a classification scheme for vertices, 
according to the configuration of polygons arranged 
around them. The abundance of different vertex types 
serves as another way to characterize disorder and iden- 
tify the manner in which polygons cluster together. In a 
perfect crystal, one would observe only vertices of type 
A. In a regular square-triangular tiling, only vertices of 
types A-D are allowed [H]. Types J and K correspond 
to the topological disclinations, while a vertex L features 
a severe pentagonal defect. To quantify the abundance of 
different vertex types, we calculate fractions Ay, defined 
as the number of vertices of a certain type Y normalized 
to the total number of particles in the polygon construc- 
tion. 

Unexpected spikes in the time dependencies of param- 
eters Ae and Ap were detected in the analysis of the re- 
cent super-heating experiment [15] , suggesting that some 
of the vertices might be metastable or exist only in a nar- 
row range of temperatures. One of the goals of our work 
is to resolve this issue. 
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Figure 1. Triangulation map of a polycrystalline 2D Yukawa 
solid, with topological defects marked as squares and trian- 
gles (a). By removing the bonds facing unusually large an- 
gles (marked by bold lines), the polygon construction (b) is 
obtained. In general there is no one-to-one correspondence 
between topological and geometrical defects. 

The polygon construction contains geometrical infor- 
mation about bond lengths and angles, as well as topo- 
logical information about the nearest-neighbour connec- 
tions. The method has also the advantage of provid- 
ing a gradation in the severity of geometrical defects. 
Quadrilaterals are the least severe, while pentagons and 
hexagons are more severe and indicate large excess vol- 
umes. Topological defects, on the other hand, provide 
only a binary measure of local orientational disorder, that 
is, at the specific location of a vertex, there either is a 
defect, or there is not. 

The gradation of defects in the polygon method allows 
for a greater sensitivity identifying and classifying disor- 
der [T5]. To characterize the state of our model system 



Figure 2. Twelve types of vertices, frequently observed in 
the polygon construction of 2D Yukawa systems. Vertices are 
classified by the number and relative distribution of polygons 
around a vertex. 

Phase transitions in two-dimensional systems are usu- 
ally identified by the sudden change in orientational or 
translational order parameters. The local orientational 
order parameter for a particle j is defined as [29 
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where 6jk is the angle between the bond connecting par- 
ticles j and k and some fixed direction. Nj is the coor- 
dination number of the particle j. The magnitude \ip6j\ 
is close to unity for a particle inside a hexagonal lattice 
but is small close to the domain walls (grain boundaries) 
or in a liquid. On the other hand, the value of a complex 
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argument arg^ej) represents the angular orientation of 
a neighbourhood or entire domain. 

The parameter ip 6 — \ (4>6j) | , that is, a magnitude of the 
averaged complex orientational order parameter, defines 
the overall orientational order of the system. In poly- 
crystalline solids, however, complex numbers ipej corre- 
sponding to the particles from different domains, tend to 
cancel in the averaging process. Therefore, approaches 
very small values in the limit of an infinite sample size. 
Another parameter can be used in such cases, namely 
tpiQi = (\ipsj\), which represents the average local orien- 
tational order of the whole system [5]. 



III. RESULTS 

We start our simulations with a defect-free lattice in 
a strongly-coupled state and the temperature of kT = 
10~ 7 . The system is then slowly heated over the pe- 
riod of time of At — 240 000, until the temperature of 
kT = 0.005 is reached. The stage of steady cooling then 
follows, restoring the temperature to its initial value (see, 
for example, Figure [8]-e). The chosen rate of heating 
is low enough to reach the equilibrium at each step of 
the simulation outside the region of a fast order-disorder 
transition. Therefore, lower rates of heating and cooling 
would not substantially change the qualitative results of 
our numerical experiments, except in the close vicinity of 
the transition, which we discuss later. 

The peak temperature is chosen to be well above the 
melting point, as found by the previous studies [TT1 150] , 
The actual kinetic temperature of the system is calcu- 
lated from observed particle velocities, kT = (vf}/2 and 
is found to fluctuate around the prescribed values of 
fcT re f, with the fluctuations being proportional to the 
temperature. We keep track of the order parameters and 
defect fractions throughout the whole cycle. In this sec- 
tion we first investigate the orientational order of the 2D 
system, and then turn to the defects and polygons. 

The initial values of both orientational order parame- 
ters ijjQ and ^lei are very close to unity, as Figure [3] shows, 
and correspond to the defect-free hexagonal lattice. The 
system exhibits a sudden loss of orientational order in 
the temperature range of kT = 0.0025-0.0027 (the cor- 
responding coupling parameter values are T = 370-400), 
which is the signature of a melting phase transition. Our 
observations are in a fair agreement with 130] , where 
phase transitions in two similar 2D systems were observed 
near the values of T = 384 and T = 415. At the end of the 
heating phase tpe drops below 0.1 and ?A|6| to the value 
of approximately 0.5. 

As the temperature is gradually lowered, a 2D Yukawa 
liquid freezes back to the hexagonal lattice. However, the 
evolution of the order parameters does not follow exactly 
the same path as in the case of melting and hysteresis 
is observed (Figure [3]) . ^6 stays low until the tempera- 
ture of kT = 0.0024 is reached, which is lower than the 
melting point. Changes in i/>| 6 | also occur at somewhat 




Figure 3. Averaged orientational order parameters ip6 an d 
■0|6| °f a 2D Yukawa system during the heating and cooling 
simulation cycle. 

lower temperatures and are not as abrupt as in the case of 
melting. As we show later, the effect of hysteresis is most 
likely a result of the finite rate of heating and cooling. 

The concentration of topological defects during the 
stage of slow heating exhibits a similar trend, with a weak 
temperature dependence before and after the transition, 
see Figure [4] A sudden proliferation of defects is ob- 
served in the range of kT 0.0025-0.0027. Changes in 
the defect concentration during the gradual cooling are 
not as abrupt and occur at somewhat lower temperatures, 
kT 0.0024-0.0022. 




Figure 4. Topological defect fraction DF as a function of the 
kinetic temperature kT. 

Four typical snapshots of particle positions during the 
transition are shown in Figure [5] Here triangles corre- 
spond to the particles with only five nearest neighbours, 
while squares mark positions of vertices with seven tri- 
angulation bonds. In a solid phase (Figure [5]-a), defects 
mostly appear as quartets, composed of two disclinations 
with five bonds and two vertices with seven neighbours. 
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Alternatively, this four-defect complex can be treated as 
a bound pair of two dislocations. Apparently, this kind of 
topological fault does not significantly change either po- 
sitional or orientational order. Larger defect complexes, 
consisting of more than four defective vertices emerge be- 
fore the melting transition. 

As it is illustrated in the second panel of the figure, 
during the melting transition (e.g. kT = 0.0026) defect 
complexes grow and spread. However, as it can be seen 
in Figure [5}b, large defect-free patches still exist, sug- 
gesting, that the transition from a defect-free lattice to 
the disordered state is not homogeneous. Bound dislo- 
cations are still found in the ordered patches, however 
are seldom seen. While some disclinations and disloca- 
tions are present, they are clearly not the main cause of 
the loss of order; most of these defects appear as parts 
of larger regions of condensed defect groups and chains. 
This observation supports the theory of grain-boundary 
induced melting [3J 31 (33], or possibly the coexistence 
of hexatic and liquid states, as described in 34 j. Fur- 
ther investigations of larger systems would be needed to 
unambiguously determine the mechanism of the melting 
transition. 




Figure 5. Typical arrangements of topological defects during 
the melting transition of a 2D Yukawa system. Triangles here 
represent particles with five closest neighbours, while squares 
mark vertices with the coordination number equal to seven. 
The corresponding temperatures are kT = 0.00245 (a), 0.0026 
(b), 0.0027 (c) and 0.0045 (d). 

Finally, at the end of the transition — Figure [5}c - 
topological defects form large inter-connected complexes, 
destroying orientational order completely. Free disloca- 
tions and disclinations can still be occasionally found in 
the larger groups of defects. The distribution of defects 
becomes homogeneous in a high-temperature liquid (Fig- 
ure [5}d) . 



As the system is slowly cooled, defect complexes and 
large defect-free regions can still be found at the tem- 
peratures as low as kT = 0.0022, together with some 
free dislocations. At lower temperatures, however, these 
complexes tend to shrink and rearrange, eventually leav- 
ing only bound and free dislocations as well as intersti- 
tial particles. The final configuration corresponds to the 
nearly perfect triangular lattice with a few free disloca- 
tions, leading to the low defect concentration and values 
of the orientational order parameter close to unity. 

Previous studies of similar systems [31] suggest that 
the effect of hysteresis might be caused by the finite 
rate of heating and cooling. We test this hypothesis by 
analysing the evolution of topological defect fraction at 
the constant prescribed temperature of fcT rc f = 0.00243, 
starting from either liquid or solid initial state. Accord- 
ing to Figure [6j at this temperature the system slowly 
switches between high and low values of the order pa- 
rameter, corresponding to the defect configurations de- 
picted in parts (a) and (c) of Figure [5j Therefore, there 
is a range of temperatures, in which ordered and dis- 
ordered states are unstable and have approximately the 
same probabilities to be observed, as reported in [3TJ [35] , 
and no hysteresis should occur in the limit of infinite sim- 
ulation time. 
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Figure 6. Time series for the topological defect fraction DF 
at the prescribed temperature of fcT re f = 0.00243. 

Let us now turn to the analysis of the polygon con- 
struction. The initial defect-free lattice corresponds to 
a triangular tiling. Therefore, the triangle is the only 
type of polygon present in the original configuration. As 
the system is continuously heated, quadrilaterals and oc- 
casional pentagons appear. As it can be seen in Figure 
[7J quadrilateral defects tend to cluster together, form- 
ing long chains or "ladders" in a solid phase (a) or 
larger patches of a distorted square lattice in a liquid 
(c). Defect-free zones arc still present during the initial 
stage of solid-liquid transition, e.g. at the temperature of 
kT = 0.0026 as depicted in Figure^-b. Solitary hexagons 
appear much later and are most abundant in the high- 
temperature Yukawa liquid (d). 
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Figure 7. Typical arrangements of geometrical defects during 
the gradual heating stage in a 2D Yukawa system: solid phase 
at t — 45 002 (a), configurations during (b) and right after (c) 
the solid-liquid transition and high-temperature liquid phase 
at t — 120 000 (d). The corresponding temperatures are kT = 
0.0018 (a), 0.0026 (b), 0.0027 (c) and 0.0045 (d). 
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The quadrilateral is a first type of geometrical defect 
to appear in a low-temperature solid, first seen near the 
temperature of kT = 5 • 10~ 4 . This is as expected, since 
the quadrilateral is the least severe geometrical defect. 
Also it is the most abundant type of defect in both solid 
and liquid states. The parameter P4, defined as a num- 
ber of quadrilateral defects normalized to 2N, increases 
steadily as the temperature rises during early stages of 
heating. As it is demonstrated in the second panel of 
Figure |8j the proliferation rate gets significantly higher 
as the temperature reaches the value of kT = 0.0025 and 
the melting transition begins. The rapid transition ends 
at around kT = 0.0027, where the order parameter fluc- 
tuates around the value of P 4 = 0.20. 

The abundance of quadrilateral defects in the two- 
dimensional Yukawa liquid does not change significantly 
as the temperature is further increased. Therefore, we 
suggest, that a 2D liquid right after the phase transition 
can be characterized by the temperature-independent 
value of the quadrilateral order parameter close to P4 — 
0.20 ± 0.01. These observations are illustrated in Figure 
[8] as time series for the order parameter (b) and temper- 
ature kT (e). 

Occasional pentagonal defects are first spotted at the 
temperature of nearly kT = 0.001. The most significant 
proliferation of these defects takes place in the range of 
kT = 0.0025-0.0027, where the value of a pentagonal 
order parameter changes from P 5 — 0.01 to P 5 = 0.03. 
After the melting transition, P5 increases almost linearly 



Figure 8. Time series for order parameters P3 (a), P4 (b), 
P5 (c) and P§ (d), characterizing the abundances of triangles, 
quadrilaterals, pentagons and hexagons in the polygon con- 
struction of the 2D Yukawa system. Time series of the kinetic 
temperature kT is presented in the last panel (e). Two verti- 
cal lines mark the beginning and ending points of the melting 
transition. 



with the temperature, at a constant rate. We may con- 
clude, that in the context of pentagonal defects, the liq- 
uid state can be characterized by values of the order pa- 
rameter P5 > 0.03 and a steady growth in a number of 
pentagons. 

In our simulations we observe only a relatively small 
number of hexagonal defects, with the highest value of an 
order parameter close to Pq « 0.012. Although the time 
and temperature dependencies of Pq are rather noisy (see 
Figure [8]-d) , some general observations can still be made. 
The most noticeable spread of hexagons starts at about 
kT = 0.0025, which roughly coincides with the start of 
the melting transition. Afterwards, Pq seems to increase 
steadily. 

Geometrical defects arrange themselves around the 
particles in a variety of ways. Some of the most frequent 
are classified in Figure [2] as vertex types A to L [T^l [H] . 
In a perfect triangular lattice, one would observe only 
the type A, where six triangles join forming a hexagon. 
In a liquid, where quadrilaterals tend to form intercon- 
nected complexes and "ladders" (Figure [7]), the number 
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of vertices B, D, G and H is expected to increase. 

As we can see, the arrangement of polygons with re- 
spect to a certain vertex can be used as an indicator 
of disorder throughout melting or freezing transitions. 
Therefore, we further investigate the evolution of vertex 
fractions Xy = Ny/N, defined as the ratio of a num- 
ber of vertices for a certain vertex type Ny to the total 
number of particles N. 

Vertices E and F are the first to appear when the 
defect-free system is slowly heated. Just as the quadri- 
lateral defects, these vertices are first observed at the 
temperature of kT w 5 • 10~ 4 . Both vertex types E and 
F contain a single quadrilateral and four or five triangles, 
and are created by removing an inner (E) or outer (F) 
bond from a vertex type A. Therefore, the abundances 
of vertices E and F have essentially identical time and 
temperature dependencies (third panel of Figure [9]). 

As the temperature rises, the order parameter Ae 
grows until the critical value of kT rj 0.0025 and the 
fraction of Ag w 0.175 is reached. As a matter of fact, 
it is approximately the same temperature that marks the 
start of the rapid melting transition, i.e. the sudden loss 
of orientational order, the rapid growth in a number of 
topological defects, squares and pentagons. As geometri- 
cal defects proliferate further at higher temperatures, the 
number of vertices E and F monotonically diminishes. 

Vertex types B and C both contain two quadrilaterals 
and three triangles, but differ in the order they are ar- 
ranged around a vertex. Two quadrilaterals of a vertex 
type B share a common edge (Figure [2]), while in a vertex 
C they are separated by a triangle or two. Vertex concen- 
trations Xb and Xq depart from zero significantly around 
kT = 0.001 and grow further with the kinetic tempera- 
ture. The order parameter Ab reaches its highest value of 
0.14 at the temperature of kT sa 0.0027, which coincides 
with the final transition to the liquid phase (Figure [9]). 
Approximately the same is true for the vertex fraction of 
the type C, which reaches its highest value Ac ~ 0.12 
close to the temperature of kT w 0.0027. As the 2D 
liquid is heated further, the concentrations of vertices B 
and C decreases. 

Vertex fractions for types D, G, H, L and I all share 
similar time and temperature dependencies (fourth panel 
of Figure [9]). They all feature a steady growth before the 
temperature of kT = 0.0025 is reached, rapid change 
throughout the melting transition and a weak tempera- 
ture dependence in a liquid state. Right after the tran- 
sition, that is, temperatures above kT = 0.0027, the 
vertex fraction D stays close to a value of Ad = 0.017, 
while those for other types fluctuate around Xq — 0.060, 
A H = 0.035, X t = 0.030, A L = 0.039. 

A few interesting results were obtained in the recent 
work [TS] , where the evolution of geometrical defects and 
vertex fractions in a laser-heated complex plasma was 
studied. It should be noted, though, that the goal of 
the study was to investigate solid super-heating and not 
the temperature dependencies of order parameters. The 
authors did not have a chance to vary the temperature 




Figure 9. Time series for vertex fractions of vertex types B-L 
and the temperature kT. The time dependence of a vertex 
fraction Xy is virtually identical to the one for vertex E and 
therefore is not shown here. Vertical lines mark the beginning 
and ending points of the rapid order-disorder transition. 

of a system gradually and a heating source was turned 
on and off abruptly. 

First, quadrilateral and pentagonal order parameters 
in a liquid state of the experimental system were found 
to be P4 « 0.19 and P5 sa 0.05. These values fully agree 
with and support the results of our simulations. Sec- 
ondly, the authors of [IS] observed sudden spikes in the 
abundances of vertex types F and E at the times when the 
heating source was abruptly turned on and off. Whether 
it is a signature of vertex metastability or a specific tem- 
perature dependence remained unclear. In the view of 
our simulation results and, specifically Figure [9j we con- 
clude that the reason behind the spikes is an intrinsic 
temperature dependence of vertex fractions B, C, E and 
F, with spikes corresponding to the temperatures lower 
than the final kinetic temperature of a liquid. 

Just as in the case of the orientational parameters 
and topological defect fractions, an evolution of order 
parameters P and vertex fractions X during the slow 
crystallization does not exactly duplicate the behaviour 
throughout the heating phase. For example, a peak in the 
time dependence of the vertex fraction E corresponds to 
the temperature of kT « 0.0022, as opposed to the tern- 



perature of kT w 0.0025 in the case melting. Critical 
values of Xb and Xq are also slightly shifted to lower 
temperatures. In a final configuration after the crystal- 
lization, there are seven quadrilateral defects and a cor- 
responding small fraction of vertices B and E (Figure [9]). 

We have repeated our simulations starting with poly- 
crystallinc configurations, obtained from a rapidly cooled 
2D Yukawa liquid. The evolution of order parameters 
throughout the initial heating phase turned out to be 
slightly configuration-dependent. For example, in a few 
cases orientational order parameters increased as the 
polycrystalline solid was heated, while at the same time 
topological defects diminished. This can be explained by 
the melting of grain boundaries and simultaneous merg- 
ing of crystalline domains. Nevertheless, the point of a 
final transition to the liquid phase was found to be the 
same as in the case of the hexagonal initial configura- 
tion, that is, kT w 0.0027. Therefore, most of the results 
presented here would also hold for the melting of a poly- 
crystalline solid. 



IV. SUMMARY 

In this contribution we report the results of Langevin 
dynamics simulations, performed to investigate two- 
dimensional melting and crystallization transitions in 
many particle Yukawa systems, such as those found in 
complex plasma experiments. To characterize the state 
of a system, we use a polygon construction method, in 
which unusually long bonds are deleted from the triangu- 
lation map of particle positions. Geometrical defects are 
identified as non-triangular polygons, while vertices are 
classified according to the type and order of polygons sur- 



rounding them. We also make use of a topological defect 
fraction and orientational order parameters as conven- 
tional tools to analyse phase transitions. 

In our simulations of the system with the constant 
value of screening strength (k = 2), the solid-to- liquid 
phase transition takes place in the coupling parameter 
range of T s=s 400-370, where rapid changes in order 
parameters and defect fractions are detected. The ori- 
entational and translational order is destroyed by the 
non-homogeneous growth of large defect complexes and 
chains, contrary to the KTHNY theory of two dimen- 
sional melting, which predicts the emergence of isolated 
dislocations and disclinations. 

To quantify the disorder, we use polygonal order pa- 
rameters P, that is, the normalized number of geomet- 
rical defects of a certain kind. It turns out, that a 
liquid phase can be characterized by the temperature- 
independent value of a quadrilateral order parameter of 
P4 w 0.20. In the context of pentagonal defects, the liq- 
uid state is characterized by the value of P5 > 0.03. 

Concentrations of vertices containing three or more 
quadrilaterals (D, G, H and I) show only a weak depen- 
dence on the temperature in the liquid state, showing the 
tendency for quadrilaterals to cluster together. Temper- 
ature dependencies of vertex fractions of the types B, C, 
E and F all feature well expressed peaks at the beginning 
(E, F) or final stages (B, C) of the solid- liquid transition 
with a very similar behaviour during the recrystallization. 
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